notebook09 - Statistical Analysis

rm(list=ls(all=TRUE))
load("~/Post-doc/Data/Total Merged Data File (July 11 2022).RData")

Data for mean velocity, angular velocity, turn and heading angles each will be analyzed by a mixed model ANCOVA with chemical concentration, and light as fixed (treatment) effects, flow velocity as the co-variate and replicate as a random (blocking affect) with individual paths as the unit of analysis.

```r
```r
library(WRS2)

Error in is.nan(CC.TotalData) : 
  default method not implemented for type 'list'

```r

legend(\topright\,                                  # Add legend to density
       legend = c(\Flow rate  0\, \Flow rate 0.6\),#, \Flow rate 3\, \Flow rate 5.9\, \Flow rate 8.9\),
       col = c(\blue\, \light blue\),#, \yellow\, \orange\, \red\),
       lty = 1)

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiRXJyb3IgaW4gc3Ryd2lkdGgobGVnZW5kLCB1bml0cyA9IFxcdXNlclxcLCBjZXggPSBjZXgsIGZvbnQgPSB0ZXh0LmZvbnQpIDogXG4gIHBsb3QubmV3IGhhcyBub3QgYmVlbiBjYWxsZWQgeWV0XG4ifQ== -->

Error in strwidth(legend, units = , cex = cex, font = text.font) : plot.new has not been called yet

<!-- rnb-output-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->



#################################################################################


                 Kernel Density Plots




##################################################################################

Turn Angles first...


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxucm0oa2QxLCBrZDIsIGtkMywga2Q0LCBrZDUsIGtkNiwga2Q3LCBrZDgpXG5cbiMjIyBLZXJuZWwgZGVuc2l0eSBwbG90cyAjIyMjIyMjIyMjIyMjIyMjIyMjIFR1cm4gQW5nbGVzXG5wYXIobWFyPSBjKDQsNCwyLDEpKVxucGFyKG1mcm93PWMoMiwyKSlcblxuIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjIyMgRmxvdyBSYXRlICMjIyMjIyMjIyMjIyMjIyMjIyMjIyMjI1xubGV2ZWxzKENDLlRvdGFsRGF0YSRGbG93LnJhdGUpXG5gYGAifQ== -->

```r
rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8)

### Kernel density plots ################### Turn Angles
par(mar= c(4,4,2,1))
par(mfrow=c(2,2))

################################ Flow Rate ########################
levels(CC.TotalData$Flow.rate)

[1] "0"   "0.6" "3"   "5.9" "8.9"

kd1 <- density(l1$turn.angle[l1$Flow.rate=="0"])
plot(kd1, xlim = c(5,175), ylim = c(0, 0.030), col='blue', lwd=2, main = "Turn Angles Vs Flow Rate, (all data no light var), Lights On")

kd2 <- density(l1$turn.angle[l1$Flow.rate=="0.6"])
lines(kd2, col='light blue', lwd=2)

kd3 <- density(l1$turn.angle[l1$Flow.rate=="3"])
lines(kd3, col='yellow', lwd=2)

kd4 <- density(l1$turn.angle[l1$Flow.rate=="5.9"])
lines(kd4, col='orange', lwd=2)

kd5 <- density(l1$turn.angle[l1$Flow.rate=="8.9"])
lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)

################################ Chlorophyll ########################
levels(CC.TotalData$Chlorophyll)

[1] "0"    "4.3"  "4.6"  "5.5"  "6.1"  "7.6"  "13.5" "19"  

kd1 <- density(l1$turn.angle[l1$Chlorophyll=="0"])
plot(kd1, xlim = c(5,175), ylim = c(0, 0.030), col='blue', lwd=2, main = "Turn Angles Vs Chlorophyll, Lights On")

kd2 <- density(l1$turn.angle[l1$Chlorophyll=="4.3"])
lines(kd2, col='light blue', lwd=2)

kd3 <- density(l1$turn.angle[l1$Chlorophyll=="4.6"])
lines(kd3, col='yellow', lwd=2)

kd8 <- density(l1$turn.angle[l1$Chlorophyll=="5.5"])
lines(kd8, col='orange', lwd=2)

kd4 <- density(l1$turn.angle[l1$Chlorophyll=="6.1"])
lines(kd4, col='red', lwd=2)

kd5 <- density(l1$turn.angle[l1$Chlorophyll=="7.6"])
lines(kd5, col='dark red', lwd=2)

kd6 <- density(l1$turn.angle[l1$Chlorophyll=="13.5"])
lines(kd6, col='purple', lwd=2)

kd7 <- density(l1$turn.angle[l1$Chlorophyll=="19"])
lines(kd7, col='black', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Chlorophyll  0", "Chlorophyll 4.3", "Chlorophyll 4.6", "Chlorophyll 5.5", "Chlorophyll 6.1", 
                  "Chlorophyll 7.6", "Chlorophyll 13.5", 'Chlorophyll 19'),
       col = c("blue", "light blue", "yellow", "orange", "red", "dark red", "purple", "black"),
       lty = 1)

################################ Guano ########################
levels(CC.TotalData$Guano)

[1] "Absent"  "Present"

kd1 <- density(l1$turn.angle[l1$Guano=="Present"])
plot(kd1, xlim = c(5,175), ylim = c(0, 0.030), col='blue', lwd=2, main = "Turn Angles Vs Guano, Lights On")

kd2 <- density(l1$turn.angle[l1$Guano=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Guano Present", "Guano Absent"),
       col = c("blue", "red"),
       lty = 1)

################################ Light ########################
levels(CC.TotalData$Light)

[1] "Absent"  "Present"

kd1 <- density(l1$turn.angle[l1$Light=="Present"])
plot(kd1, xlim = c(5,175), ylim = c(0, 0.030), col='blue', lwd=2, main = "Turn Angles Vs Light")

kd2 <- density(l0$turn.angle[l0$Light=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)

Now looking at velocity in relation to flow instead….

### Kernel density plots ################### Swimming velocity
par(mar= c(4,4,2,1))
par(mfrow=c(2,2))
#range(CC.TotalData$vel.flow)
#plot(CC.TotalData$vel.flow)
CC.TotalData$log.vel.flow <- log10(CC.TotalData$vel.flow)

Warning: NaNs produced

#str(CC.TotalData)
################################ Flow Rate ########################
levels(CC.TotalData$Flow.rate)

[1] "0"   "0.6" "3"   "5.9" "8.9"

kd1 <- density(log10(l1$vel.flow[l1$Flow.rate=="0"]))

Warning in density(log10(l1$vel.flow[l1$Flow.rate == "0"])) :
  NaNs produced
Error in density.default(log10(l1$vel.flow[l1$Flow.rate == "0"])) : 
  'x' contains missing values

Now for headings…..

rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8)

Warning in rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8) :
  object 'kd2' not found
Warning in rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8) :
  object 'kd3' not found
Warning in rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8) :
  object 'kd4' not found
Warning in rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8) :
  object 'kd5' not found
Warning in rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8) :
  object 'kd7' not found
Warning in rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8) :
  object 'kd8' not found

Now for pitch….

### Kernel density plots ################### Pitch
par(mar= c(4,4,2,1))
par(mfrow=c(2,2))

CC.TotalData$pitch.perfect <- CC.TotalData$pitch/(2*pi)*360
################################ Flow Rate ########################
levels(CC.TotalData$Flow.rate)

[1] "0"   "0.6" "3"   "5.9" "8.9"

kd1 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Flow.rate=="0"])
plot(kd1, xlim = c(-90,90), ylim = c(0, 0.06), col='blue', lwd=2, main = "Pitch Vs Flow Rate")

kd2 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Flow.rate=="0.6"])
lines(kd2, col='light blue', lwd=2)

kd3 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Flow.rate=="3"])
lines(kd3, col='yellow', lwd=2)

kd4 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Flow.rate=="5.9"])
lines(kd4, col='orange', lwd=2)

kd5 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Flow.rate=="8.9"])
lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)

################################ Chlorophyll ########################
levels(CC.TotalData$Chlorophyll)

[1] "0"    "4.3"  "4.6"  "5.5"  "6.1"  "7.6"  "13.5" "19"  

kd1 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Chlorophyll=="0"])
plot(kd1, xlim = c(-90,90), ylim = c(0, 0.06), col='blue', lwd=2, main = "Pitch Vs Chlorophyll")

kd2 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Chlorophyll=="4.3"])
lines(kd2, col='light blue', lwd=2)

kd3 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Chlorophyll=="4.6"])
lines(kd3, col='yellow', lwd=2)

kd8 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Chlorophyll=="5.5"])
lines(kd8, col='orange', lwd=2)

kd4 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Chlorophyll=="6.1"])
lines(kd4, col='red', lwd=2)

kd5 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Chlorophyll=="7.6"])
lines(kd5, col='dark red', lwd=2)

kd6 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Chlorophyll=="13.5"])
lines(kd6, col='purple', lwd=2)

kd7 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Chlorophyll=="19"])
lines(kd7, col='black', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Chlorophyll  0", "Chlorophyll 4.3", "Chlorophyll 4.6", "Chlorophyll 5.5", "Chlorophyll 6.1", "Chlorophyll 7.6", "Chlorophyll 13.5", 'Chlorophyll 19'),
       col = c("blue", "light blue", "yellow", "orange", "red", "dark red", "purple", "black"),
       lty = 1)

################################ Guano ########################
levels(CC.TotalData$Guano)

[1] "Absent"  "Present"

kd1 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Guano=="Present"])
plot(kd1, xlim = c(-90,90), ylim = c(0, 0.06), col='blue', lwd=2, main = "Pitch Vs Guano")

kd2 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Guano=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Guano Present", "Guano Absent"),
       col = c("blue", "red"),
       lty = 1)

################################ Light ########################
levels(CC.TotalData$Light)

[1] "Absent"  "Present"

kd1 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Light=="Present"])
plot(kd1, xlim = c(-90,90), ylim = c(0, 0.06), col='blue', lwd=2, main = "Pitch Vs Light")

kd2 <- density(CC.TotalData$pitch.perfect[CC.TotalData$Light=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)

For each level of Flow

For each level of Chlorophyll

### Kernel density plots ################### Turn Angles
par(mar= c(4,4,2,1))
par(mfrow=c(3,3))

################################ Flow Rate ########################
levels(CC.TotalData$Flow.rate)

[1] "0"   "0.6" "3"   "5.9" "8.9"

kd1 <- density(ChlA0$turn.angle[ChlA0$Flow.rate=="0"])
plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 0")

kd2 <- density(ChlA0$turn.angle[ChlA0$Flow.rate=="0.6"])
lines(kd2, col='light blue', lwd=2)

kd3 <- density(ChlA0$turn.angle[ChlA0$Flow.rate=="3"])
lines(kd3, col='yellow', lwd=2)

kd4 <- density(ChlA0$turn.angle[ChlA0$Flow.rate=="5.9"])
lines(kd4, col='orange', lwd=2)

kd5 <- density(ChlA0$turn.angle[ChlA0$Flow.rate=="8.9"])
lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)

##kd1 <- density(ChlA4.3$turn.angle[ChlA4.3$Flow.rate=="0"])
##plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 4.3")

kd2 <- density(ChlA4.3$turn.angle[ChlA4.3$Flow.rate=="0.6"])
plot(kd2, xlim = c(5,175), ylim = c(0, 0.02), col='light blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 4.3")
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA4.3$turn.angle[ChlA4.3$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA4.3$turn.angle[ChlA4.3$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA4.3$turn.angle[ChlA4.3$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)

kd1 <- density(ChlA4.6$turn.angle[ChlA4.6$Flow.rate=="0"])
plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 4.6")

kd2 <- density(ChlA4.6$turn.angle[ChlA4.6$Flow.rate=="0.6"])
##plot(kd2, xlim = c(5,175), ylim = c(0, 0.02), col='light blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 4.6")
lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA4.6$turn.angle[ChlA4.6$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA4.6$turn.angle[ChlA4.6$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA4.6$turn.angle[ChlA4.6$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(ChlA5.5$turn.angle[ChlA5.5$Flow.rate=="0"])
plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 5.5")

##kd2 <- density(ChlA5.5$turn.angle[ChlA5.5$Flow.rate=="0.6"])
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA5.5$turn.angle[ChlA5.5$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA5.5$turn.angle[ChlA5.5$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA5.5$turn.angle[ChlA5.5$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



##kd1 <- density(ChlA6$turn.angle[ChlA6$Flow.rate=="0"])
##plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 6")

kd2 <- density(ChlA6$turn.angle[ChlA6$Flow.rate=="0.6"])
plot(kd2, xlim = c(5,175), ylim = c(0, 0.02), col='light blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 6")

##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA6$turn.angle[ChlA6$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA6$turn.angle[ChlA6$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA6$turn.angle[ChlA6$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(ChlA7$turn.angle[ChlA7$Flow.rate=="0"])
plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 7")

##kd2 <- density(ChlA7$turn.angle[ChlA7$Flow.rate=="0.6"])
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA7$turn.angle[ChlA7$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA7$turn.angle[ChlA7$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA7$turn.angle[ChlA7$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



##kd1 <- density(ChlA13$turn.angle[ChlA13$Flow.rate=="0"])
##plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 13")

kd2 <- density(ChlA13$turn.angle[ChlA13$Flow.rate=="0.6"])
plot(kd2, xlim = c(5,175), ylim = c(0, 0.02), col='light blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 13")

##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA13$turn.angle[ChlA13$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA13$turn.angle[ChlA13$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA13$turn.angle[ChlA13$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)



kd1 <- density(ChlA19$turn.angle[ChlA19$Flow.rate=="0"])
plot(kd1, xlim = c(5,175), ylim = c(0, 0.02), col='blue', lwd=2, main = "Turn Angle Vs Flow Rate, Chlorophyll = 19")

##kd2 <- density(ChlA19$turn.angle[ChlA19$Flow.rate=="0.6"])
##lines(kd2, col='light blue', lwd=2)

##kd3 <- density(ChlA19$turn.angle[ChlA19$Flow.rate=="3"])
##lines(kd3, col='yellow', lwd=2)

##kd4 <- density(ChlA19$turn.angle[ChlA19$Flow.rate=="5.9"])
##lines(kd4, col='orange', lwd=2)

##kd5 <- density(ChlA19$turn.angle[ChlA19$Flow.rate=="8.9"])
##lines(kd5, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Flow rate  0", "Flow rate 0.6", "Flow rate 3", "Flow rate 5.9", "Flow rate 8.9"),
       col = c("blue", "light blue", "yellow", "orange", "red"),
       lty = 1)


###########################################################################

#############################  Light  ######################################

###########################################################################

par(mar= c(4,4,2,1))

par(mfrow=c(3,3))

kd1 <- density(ChlA0$turn.angle[ChlA0$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 0")

kd2 <- density(ChlA0$turn.angle[ChlA0$Light=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)

kd1 <- density(ChlA4.3$turn.angle[ChlA4.3$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 4.3")

##kd2 <- density(ChlA4.3$turn.angle[ChlA4.3$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)

kd1 <- density(ChlA4.6$turn.angle[ChlA4.6$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 4.6")

##kd2 <- density(ChlA4.6$turn.angle[ChlA4.6$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)

kd1 <- density(ChlA5.5$turn.angle[ChlA5.5$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 5.5")

##kd2 <- density(ChlA5.5$turn.angle[ChlA5.5$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)

kd1 <- density(ChlA6$turn.angle[ChlA6$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 6")

##kd2 <- density(ChlA6$turn.angle[ChlA6$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)

kd1 <- density(ChlA7$turn.angle[ChlA7$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 7")

kd2 <- density(ChlA7$turn.angle[ChlA7$Light=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)



kd1 <- density(ChlA13$turn.angle[ChlA13$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 13")

##kd2 <- density(ChlA13$turn.angle[ChlA13$Light=="Absent"])
##lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)


kd1 <- density(ChlA19$turn.angle[ChlA19$Light=="Present"])
plot(kd1,xlim = c(5,175), ylim = c(0, 0.025), col='blue', lwd=2, main = "Turn Angles Vs Light, Chlorophyll = 19")

kd2 <- density(ChlA19$turn.angle[ChlA19$Light=="Absent"])
lines(kd2, col='red', lwd=2)

legend("topright",                                  # Add legend to density
       legend = c("Light Present", "Light Absent"),
       col = c("blue", "red"),
       lty = 1)

For each level of Flow (now for velocity)

### Kernel density plots ################### Velocity
par(mar= c(4,4,2,1))
par(mfrow=c(3,2))

################################ Chlorophyll ########################

kd1 <- density(log10(f0$vel.flow[f0$Chlorophyll=="0"]))

Warning in density(log10(f0$vel.flow[f0$Chlorophyll == "0"])) :
  NaNs produced
Error in density.default(log10(f0$vel.flow[f0$Chlorophyll == "0"])) : 
  'x' contains missing values

For each level of Chlorophyll (velocity)

### Kernel density plots ################### velocity
par(mar= c(4,4,2,1))
par(mfrow=c(3,3))

################################ Flow Rate ########################
kd1 <- density(log10(ChlA0$vel.flow[ChlA0$Flow.rate=="0"]))

Warning in density(log10(ChlA0$vel.flow[ChlA0$Flow.rate == "0"])) :
  NaNs produced
Error in density.default(log10(ChlA0$vel.flow[ChlA0$Flow.rate == "0"])) : 
  'x' contains missing values

Now for each level of Flow (for headings)

median(f0$heading.pi[f0$Chlorophyll=="0"])

[1] 0

Now for each level of Chlorophyll (for headings)